home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / simpson.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  173 lines
  1. ; $Id: simpson.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ; Copyright (c) 1991-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. function simprec1,funct,a, b,y, level, tol, count, S0, p, complex
  7. return,0
  8. end
  9.  
  10.  
  11. function simprec1, funct, a, b, y, level, tol, count, S0, p, complex
  12. ON_ERROR,2
  13. common simp2, reclev
  14. maxlev = 10
  15. if level > maxlev THEN BEGIN
  16. reclev=1
  17.  count = 0
  18.  m = (a+b)/2
  19.  if p THEN BEGIN
  20.    v = CALL_Function(funct,m)
  21.    if complex eq 0 THEN $
  22.       plots, [m, v], psym = 2  $
  23.    else BEGIN
  24.       plots, [m, float(v)], psym = 2
  25.       plots, [m, imaginary(v)], psym=4
  26.    ENDELSE
  27.    wait,1
  28.  ENDIF
  29.  
  30.  return, S0
  31. ENDIF ELSE BEGIN
  32.  
  33. h = (b - a)/2.
  34. m = (a + b)/2.
  35. int =  [a + h/2., b - h/2.]
  36. if complex eq 0 then y1 = fltarr(2) else y1=complexarr(2)
  37. for i = 0,1 do y1(i) = CALL_FUNCTION(funct, int(i))
  38. count = 2
  39.  
  40. if p ne 0 THEN BEGIN
  41.   if complex eq 0 THEN $
  42.     plots, int, y1, psym=2 $
  43.   else BEGIN
  44.     plots, int, float(y1), psym = 2
  45.     plots, int, imaginary(y1), psym=4
  46.   ENDELSE
  47.   wait, 1
  48. ENDIF
  49.   
  50.   
  51. S1 =  h*(y(0) + 4*y1(0) + y(1))/6.
  52. S2 =  h*(y(1) + 4 * y1(1) + y(2))/6.
  53. S = S1 + S2
  54.  
  55. if abs(S - S0) gt tol * abs(S) THEN BEGIN
  56.    S1 = simprec1(funct,a,m, [y(0), y1(0), y(1)], level + 1, tol/2, c1, S1, $
  57.                  p, complex)
  58.    S2 = simprec1( funct,m, b, [y(1), y1(1), y(2)], level + 1, tol/2, c2, S2, $
  59.                   p, complex)
  60.    count = count + c1 + c2
  61.    return, S1 + S2
  62. ENDIF ELSE return, S
  63.  
  64. ENDELSE
  65. END
  66.  
  67.  
  68.  
  69. function simpson, funct, a, b, count,  tol=tol, plot_it = plot_it, $
  70.                    complex = complex
  71. ;+
  72. ; NAME:
  73. ;    SIMPSON
  74. ;
  75. ; PURPOSE:
  76. ;    Numerically approximate the definite integral of a function 
  77. ;    with limits [A, B].
  78. ;
  79. ; CATEGORY:
  80. ;    Mathematical functions, general.
  81. ;
  82. ; CALLING_SEQUENCE:
  83. ;    Result = SIMPSON(Funct, A, B, Count, TOL = Tol)
  84. ;
  85. ; INPUTS: 
  86. ;    Funct:    A character string that contains the name of the function
  87. ;        to be integrated.  The user should write this function
  88. ;        to accept a single scalar argument.
  89. ;
  90. ;    A:    The lower limit of the integral.
  91. ;
  92. ;    B:    The upper limit of the integral.
  93. ;
  94. ; KEYWORD PARAMETERS:
  95. ;    TOL:    The error tolerance. The default is .001.    
  96. ;
  97. ;     COMPLEX:    Set this keyword if Funct returns complex values.
  98. ;
  99. ;     PLOT_IT:    Set this keyword to plot the points where the integrand is 
  100. ;        evaluated.
  101. ;
  102. ; OUTPUT:
  103. ;    SIMPSON returns the approximation to the integral of Funct for the 
  104. ;    limits [A, B].
  105. ;
  106. ; OUTPUT PARAMETERS:
  107. ;    Count:    On return, this variable contains the number of function 
  108. ;        evaluations needed to approximate the integral.
  109. ;
  110. ; COMMON BLOCKS:
  111. ;    SIMP2
  112. ;
  113. ; SIDE EFFECTS:
  114. ;    None.
  115. ;
  116. ; RESTRICTIONS:
  117. ;    None.
  118. ;
  119. ; PROCEDURE:
  120. ;    This function uses the recursive adaptive Simpson's rule.
  121. ;
  122. ; MODIFICATION HISTORY:
  123. ;    1991, Ann Bateson.
  124. ;
  125. ;-
  126.  
  127.  
  128. On_Error,2
  129.  
  130. common simp2, reclev
  131.  
  132. reclev = 0
  133. if KEYWORD_SET(tol) eq 0 THEN tol = .001
  134. if KEYWORD_SET(complex) eq 0 THEN complex = 0
  135. if a eq b THEN return,0
  136. if KEYWORD_SET(plot_it) eq 0 THEN p = 0 $
  137. else BEGIN
  138.   p = 1
  139.   m = (b+a)/2
  140.   x = [m, b, a + findgen(11) * (b-a)/10.]
  141.   if complex eq 0 THEN v = fltarr(13) else v = complexarr(13)
  142.  
  143.   for i = 0,12 do v(i) = CALL_FUNCTION(funct,x(i))
  144.  
  145.   if complex eq 0 THEN BEGIN
  146.      plot, x, v, /nodata  
  147.      plots, x(0:2), v(0:2), psym=2
  148.   ENDIF ELSE BEGIN
  149.      v1 = imaginary(v)
  150.      v  = float(v)
  151.      plot, [x,x], [v, v1], /nodata
  152.      plots,x(0:2), v(0:2), psym=2
  153.      plots, x(0:2), v1(0:2), psym=4
  154.   
  155.   ENDELSE
  156.   wait,1
  157. ENDELSE
  158.  
  159. level = 1
  160. m = (a + b)/2
  161.  
  162. if complex eq 0 THEN y = fltarr(3) else y = complexarr(3)
  163. x = [ a, m, b]
  164. for i = 0,2 do  y(i) = CALL_FUNCTION(funct, x(i))
  165.  
  166. int = simprec1(funct, a, b, y, level, tol, count, 1.e+30, p, complex)
  167.  
  168. if reclev ne 0 THEN $
  169.   print, "Recursion limit has been exceeded. Beware singularity"
  170. count = count + 3
  171. return, int
  172. END
  173.